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ABSTRACT 

StarBench is a project focused on benchmarking and validating different star-formation 
and stellar feedback codes. In this first StarBench paper we perform a comparison study of 
the D-type expansion of an H II region. The aim of this work is to understand the differences 
observed between the twelve participating numerical codes against the various analytical ex¬ 
pressions examining the D-type phase of H II region expansion. To do this, we propose two 
well-defined tests which are tackled by ID and 3D grid- and SPH- based codes. The first test 
examines the ‘early phase’ D-type scenario during which the mechanical pressure driving the 
expansion is significantly larger than the thermal pressure of the neutral medium. The second 
test examines the Tate phase’ D-type scenario during which the system relaxes to pressure 
equilibrium with the external medium. Although they are mutually in excellent agreement, all 
twelve participating codes follow a modified expansion law that deviates significantly from 
the classical Spitzer solution in both scenarios. We present a semi-empirical formula com¬ 
bining the two different solutions appropriate to both early and late phases that agrees with 
high-resolution simulations to < 2%. This formula provides a much better benchmark solu¬ 
tion for code validation than the Spitzer solution. The present comparison has validated the 
participating codes and through this project we provide a dataset for calibrating the treatment 
of ionizing radiation hydrodynamics codes. 

Key words: (ISM:) H II region - ISM: kinematics and dynamics - methods: numerical - 
hydrodynamics - ISM: bubbles - galaxies: ISM 


1 INTRODUCTION 
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to complex systems beyond the capability of analytic calculation, 
and also provide validation for the physical correctness of proposed 
theoretical models. However, computational models have their own 
sources of error. It is therefore important to validate their accuracy, 
where possible, if their results are to have any credence. Different 
sources of error require different validation techniques, but use of 
code inter-comparisons and resolution studies for well-defined test 
problems is particularly powerful in this regard. Previou s suc h com¬ 
pariso ns have included photoio nization modelling dPeauignot et al.l 
l200ll) . PDR code co mparison dRollig et al. 2001) . and cosmologi¬ 
cal code comparison l llliev et all200dl2009h . 

As well as the standard tests which are widely applied to com¬ 
putational hydrodynamics codes, such as Sod shoc k tubes (ISod 
1978) and the Sedov-Taylor point explosion ( Sedovl ll959l : iTavlon 


19501) . it is important to develop test problems more specific to the 


domain of study, which test broader physics couplings. The Star- 
Bench project seeks to develop a collection of well-understood 
validation cases for the modelling of star-forming regions. To date, 
two workshops have been held (Exeter, April 2013; Bonn, Septem¬ 
ber 20140, in which groups specializing in numerical star forma¬ 
tion have gathered to set the basis of standard and well-defined 
tests and explore the differences observed between their codes. The 
present paper reports the results of the first test investigated by the 
StarBench project: the D-type expansion of an H II region. 

The problem of the expansion of an H II region into a uni¬ 
form neutral medium has a long history in theoretical and numer¬ 
ical studies of the interstellar medium (ISM). The ionizing pho¬ 
tons (hv > 13.6 eV) produced by a hot, massive star ionize and 
heat a volume of gas around the star, which then begins to expand 
supersonically into the surrounding medium due to the overpres¬ 
sure. The photoionized region is bounded by an ionization front , 
which separates th e ionized from the neutral gas. Kahn! Jl954h 
and lAxfordl Jl96lh showed analytically that the ionization front 
relatively quickly switches from R-type to D-type and drives a 
strong shock into the neutral ISM. This behaviour was verified 
by early numerical invest igati ons of H II re gion expansion con¬ 
ducted by [Mathewsl (1965) and Laskerl dl966h . These early analytic 
and n umerical calculations are reviewed by iMathews & O’Delll 
( 1969 ). and subsequent developments, including the first multi¬ 
dimensional simula tions in uniform and non-uniform media, are 
reviewed bv I Yorkel d 1986 ). An analytic solution for the expansion 
rate of the H II region during the D-type phase as a function of time 
was derived by Sp it zer ( 1978:), which is presented more lucidly in 
the book of iDvson & Williams ( 1980 ). This has become the stan¬ 
dard solution with which numerical results are compared. 

Most modern code developers use the spherical expansion 
of an H II region as a standard test problem to validate their 


(e.g. 

Mellema et al. 2006a: 

iKrumholz et al. 2007: Mackev & Lim 

I 2010 

,2011 HTremblin et al.l 

l2012al: IRaga et al.ll2012al). This nrob- 


lem involves a complicated combination of fluid dynamics, radia¬ 
tive transfer, microphysical heating, cooling, ionization and recom¬ 
bination, and differences between results obtained with different 
codes and algorithms have not so far been addressed in any depth 
in the literature. There is therefore a pressing need of setting up 
benchmarking tests to examine the behaviour of the different codes 
in reproducing known anal ytical expressions. 

Although the lSpitze f 11978 ) analytical solution is widely used 


for code testing, alternative analytical solutions for the expan¬ 
sion of an H II region, applicable at different stages of the evo¬ 
lution and based on different assumptions, have been devel oped 
more recently, suc h as those of iHosokawa & Inutsukal (120061) and 
IRaga et all d2012 a.b). It is an aim of the present paper to determine 
which, if any, of the available analytical solutions is most suitable, 
and thus to provide a benchmark useful for any other radiation hy¬ 
drodynamics code. For instance, in the first attempts of the present 
test, several discrepancies between the codes have been identified 
resulting from bugs or incorrect assumptions which have all been 
fixed throughout this work. This demonstrates the value of such a 
comparison in improving confidence in the accuracy of these codes 
when applied to more complex scenarios. 

The interest of the authors of the present paper in establishing 
a reliable benchmark test for expanding H II regions stems from 
the large body of work, much of it done in the last decade, on the 
influence of expanding H II regions at and below giant molecular 
cloud (GMC) scales, intended to explore explicitly the role of pho¬ 
toionisation in the evolution and destruction/dispersal of GMCs and 
the regulation of star-formation rates and efficiencies. H II region 
feedback has been invoke d both as a trigger of star formation (e.g. 
lElmegreen & Ladal Il977l). and as the main agent re sponsible for 
terminating it (e.g. lMcKee et al.l!l984lMatznen!2002l) . Many simu¬ 
lations have sought to model one or both of these scenarios in order 
to establish whether the influence of photoionising feedback on the 
galactic star-formation process is overall positive or negative. 

For example, the gravitational instability of the dense shells 
driven by H II regions, often referr ed to as the ‘collect an d 
collapse process has been modelled b y Elmegreen & Ladal dl977l) : 
I Whit worth et al.l dl994al) : iDale et all ( 20071) . who show that, in 
principle, it should be an efficient means of triggering star for¬ 
mation in even smooth media. Other authors have instead studied 
the influence of substructure encountered by the radiation field, 


driven implosion 

model (Bertoldi 

19891: Lefloch&La 

zareffl 

1994 Mellema et al.l 19981; Kesse 

-Devnet & Burked 

20031: 

Hennev et al. |2009|; Mackev & Lim 

20111: Bisbas et al] 

20TT1; 

Haworth & Harries 

2012i). or of more generalised fractal sub- 
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l201lLl20lil20l3h . 


The negative effects of H II region expansion on clouds 
alrea dy in a dynamica l stat e have also receiv ed much atten¬ 
tion. IPeters et all d201Q[) and IPeters et all d201ll) model H II re¬ 
gions driven by roughly centrally-located massive stars on rotat¬ 
ing GMCs, finding them rather ineffective at preventing the for¬ 
mation of stars. The general problem of the interaction of ex- 


bv Mellema et al. ( 2006al): Gritschneder et al. ( 20091 

): Arthur et al. 

(2011 

): iTremblin et al. ( 

I2012bl): IDale et al. 12012 

1): IColm et al. 


( 20131) : iBoneberg et al.l d2015l) . most of whom conclude that H II 


regions are able to drive turbulence, and are destructive to GMCs 
of ~ 10 4 M©. 

This paper is organized as follows. In Section |2] we present 
a theoretical background of the D-type expansion of an H II re¬ 
gion and discuss the different analytical solutions to this problem. 
In Section [3] we present the numerical codes participating in this 
test, highlighting in brief the methods used to propagate the ioniz¬ 
ing radiation. In Section |4| we present the initial conditions of the 
1 StarBench-I: http://www.astro.ex.ac.uk/people/haworth/workshop_bench/ind^M® narking test and in Section[5]we show the results. We sum- 
StarBench-II: https://www.astro.uni-bonn.de/sb-ii/index.html marize and conclude in Section[6] 
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2 THEORETICAL BACKGROUND 


A hot, massive star emits large numbers of extreme ultraviolet 
(EUV) photon with energy ho larger than the ionization poten¬ 
tial of hydrogen (13.6eV). The ionizing photons interact with the 
surrounding neutral medium and ionize a volume of gas within 
which the rate of ionizations is nearly balanced by recombinations 
with the remaining small flux ionizing the neutral material passing 
through its surface. The kinetic energy of the ejected photoelectrons 
heats the gas within the photoionized volume by collisions but ef¬ 
ficient cooling due to collisionally excited line radiation of metal 
ions such as doubly ionized oxygen means that temperatures inside 
H II regions are approximately uniform and typically ~ 10 4 K (see 
IOsterbrocklll989L for details). In the analysis that follows, we will 
ignore the detailed physical processes and make a series of simpli¬ 
fying assumptions. 

Consider a spherically symmetric cloud of radius R c \, of to¬ 
tal mass M c i and of uniform density p c \ consisting of pure atomic 
hydrogen with a uniform temperature T c \. Let us place a radiation 
source at the centre of the cloud which defines the origin of a Carte¬ 
sian coordinate system. The source emits A/* LyC Lyman continuum 
ionizing photons per unit time and we assume that the ionizing pho¬ 
tons are monochromatic with energy ho = 13.6 eV. We assume a 
photoionization cross-section of a = 6.3 x 10“ 18 cm 2 and we use 
the recombination coefficient, a B , int o excited stages only by in¬ 
voking the on-the-spot approximation (lOsterbrockll 19891 . known as 
the ‘case-B’ recombination coefficient). Assuming an isothermal 
H II region at Ti = 10 4 K the case-B recombination coefficient is 
taken to be a B ~ 2.7 X 10 -13 cm 3 s -1 . Hereafter, the indices T 
and ‘o’ shall denote the ionized and the neutral medium, respec¬ 
tively. 

IStromgrenl il939t) was the first to show that the transition zone 
between the ionized gas and the surrounding neutral medium oc¬ 
curs over a very short distance compared to the dimensions of the 
H II region. This transition zone is the ionization front and we can 
treat it as a sharp discontinuity. The distance over which the degree 
of ionization changes from 90% to 10% is given by 


Aflst ~ ~ 1.72 x 10 -4 


p 0 cr 


pc 


Po 


IQ - 20 gem -3 


( 1 ) 


(e.g. 


Whitw orth !200(f) where m p represents the proton mass. 

Kahnl dl954l) studied in detail the propagation of an ioniza¬ 


tion front into a neutral medium when an ionizing star suddenly 
switches on. At early times most of the Lyman continuum pho¬ 
tons ionize additional gas beyond the instantaneous position of the 
ionization front. Thus the H II region expands rapidly at highly su¬ 
personic speed relative to the sound speed of the ionized gas. In 
this first phase the i oniz ation front is called R-type (R=Rarefied). 
As shown bv IStromgrenl Jl939h . the A// yC photons emitted by the 
central source will ionize a spherical region of radius 


-Rst — 


3-A/~LyC m p 
47T a B p§ 


1/3 


( 2 ) 


where p 0 = p c \ and m p is the proton mass. Hereafter, we refer to 
the spherical region of radius Rst as the initial Stromgren sphere 
and the radius .Rst as the initial Stromgren radius. The R-type 
phase of expansion terminates once the ionization front reaches this 
initial Stromgren radius. The timescale for this first phase is of or¬ 
der the recombination time 


m p 

P° 


19.6 yrs 


10 - 2 ° gem -3 


(3) 


The large temperature difference between the two regions re¬ 
sults in a large difference in thermal pressure, and the H II region 
expands. In this second phase, which is called D-type (D=Dense), 
the ionization front propagates at subsonic speed relative to the ion¬ 
ized gas but at supersonic speed relative to the neutral gas. There¬ 
fore, it is preceded by a shock front which sweeps up a dense shell 
of neutral gas. The shell is bounded on its inside by the ionization 
front and on its outside by the shock front. 

Within the ionized region we assume that at all times there is a 
balance between the ionizing photons produced by the star and the 
recombination events. Therefore, 

Kyc = ( -^ £l -<*BdV, (4) 

y J m p m e 

where p e and m e are the electron density and mass respectively. 
If we assume that i) the medium is entirely ionized, ii) p e = 
Pim e /m p , and iii) that the ionization front has a sharp edge, the 
above integral leads to 


•^LyC — 


47T 

~3~ 





which, when combined with Eqn. gives 


(5) 


pi(t) = po | 


-Rst 

R IF {t) 


3/2 


( 6 ) 


where R 1F is the extend of the ionized region. As the H II region 
expands, it ionizes more neutral gas and its mass increases in time. 
The time-dependent mass of the ionized region, Mft) is given by: 

M i (t) = ^- Po R 3 s i 2 R 3 1 l 2 (t). (7) 


The thermal pressure of the ionized gas in an H II region 
(which drives the expansion) matches approximately with the ther¬ 
mal pressure of the neutral gas in the shell between the ionization 
front and the shock front. 


2.1 Early time behaviour 

2.1.1 Spitzer approximation 


The first analytical attempt to model t he D-type expansion of an 
H II region was perfor med bvISpitzed (1 19781) wh ich was investi¬ 
ga ted in more depth b y lDvson & William j ( 198(3) . As pointed out 
by iRaga et all d2012ah . by assuming the thin shell approximation 
and by equating the pressure of the neutral gas in the shell between 
the ionization front and the shock front with the ram pressure of 
the undisturbed neutral gas as it is swept up by the shock front we 
obtain: 


1 dR Sp (t) _ J 

r Rst i 

1 3/4 Mi To 1 

r Rst \ 

ci dt ] 


1 PoTi) 

[Rs P (t)j 


- 3/4 


( 8 ) 


where p is the mean molecular weight of the gas. The above equa¬ 
tion shall be referred to as ‘Raga-T throughout the present paper. 
The ratio piT 0 /p 0 Ti found on the right-hand-side of Eqn. 
is generally small (~ 1/200). At early times the term 

777^ { R^ S (t) } can ne gl ecte d, therefore Eqn. ([8} leads to 
the so-called Spitzer solution: 


Rsp(t) — Rst (1 + 


7 at \ 
4 Rst J 


4/7 


(9) 
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2.1.2 Hosokawa-Inutsuka approximation 


2.2.2 Raga’s extension of Hosokawa-Inutsuka (Raga-II) 


A different approach describing the expan sion of an HII region was 
provided bv lHosokawa & Inutsuka 12006 ) who derived the position 
of the ionization front in time directly from the equation of motion 
of the expanding shell. Independently, iRaga et afl ( 2012b ) argued 
that the differential Eqn. © does not incorporate the inertia of the 
shocked neutral gas which is created during the expansion of the 
H II region. To include the inertia, we may write the equation of 
motion of the mass M within the shocked shell as 

f t (MRm) = 4t rR 2 m {Pi - P a ), ( 10 ) 

where P is the thermal pr essure of the gas and Rm is t he position 
of the shell as modelled bv lHosokawa & Inutsukal (120061 ‘HI’). The 
term Pi — P Q corresponds to the net thermal pressure taking into 
account the pressure acting from the neutral gas onto the ionized 
gas. Using Eqn. © and where Pif (t) represents now the position 
of the shocked shell Phi, we obtain the second-order non-linear 
differential equation 


Phi + 


(£) 


Q Z ? 3 / 2 r 2 

•2 _ on st q 
Jx hi — 


R\ 


,5/2 


3cp 

Phi 


( 11 ) 


The above equation shall be referred to as ‘Raga-II’ throughout the 
present paper. Equation cld represents the equation of motion of 
the expanding shell when its inertia is taken into account. Solving 
the above leads to 


Phi (f) 



4 

3 


gjjT 0 
2/ioP' 


( 12 ) 


At early times we may neglect the gL]T 0 / p 0 T{ term. This corre¬ 
sponds to the assumption of P Q = 0 in Eqn. dTOt . Equation IT2l has 
therefore the analytical solution 

• (i3> 

This equation was first presented bv lHosokawa & Inutsukal (2006) 
and it differs from the known Spitzer expression (Eqn. B by a 
^/4/3 factor. This factor arises from our inclusion of the inertia 
of the shocked gas due to its own movement leading to a slightly 
faster expansion than that obtained from Eqn. {9]). 


2.2 Late time behaviour 

2.2. 1 Raga’s extension of Spitzer ( Raga-I) 

At later times the term j R ^ st ^ j increases and eventu¬ 
ally the H II region stagnates at t — f STAG which is defined by 
Ps P (f STAG ) = 0. By this time the H II region is in pressure equi¬ 
librium, thus it does not expand further. The stagnation radius is 
(from Eqn[ 8 j 


/ \ 4/3 

Kstag = (^J ifet, (14) 

the density of the ionized gas is (from Eqn[ 6 ]) 



and the total ionized mass is (from Eqn0 


Mi(t STAG ) = yi&Po (A) • ( 16 ) 


At later times Phi becomes large, and therefore the two terms con¬ 
tained in the square root of Eqn. 03 become comparable. Eventu¬ 
ally at t — t STAG in which Pm(f STAG ) = 0 we obtain the stagna¬ 
tion radius 

«,™ = fls,(|) 2,, (|) 4 ' 3 , <>7) 

the density of the ionized medium 



and the total ionized mass 


4 /o\ 8/3 / \ 8/3 

M(i STAG ) = y Rlpo (§) (f) • (19) 


3 NUMERICAL CODES 

The method^ which the present numerical codes used to ac¬ 
count for hydrodynamics a r e either gr i d-based techniques (e.g . 
IColella & Woodwardl ll984l : iRoel 1 19861 : iBerger & Colellal Il989l) 
or smoothed particle hydrodynamics-ba sed techniques (SPH; e.g. 
iLucvl 1 1 97~7l : iGingold &Monaghanll 19771) . Here we present a brief 
review of the participating codes and their radiative transfer capa¬ 
bilities. All codes are summarized in Table [l] 


3.1 aquiline, glide (R. Williams) 

AQUILINE is a simplified version of AQUALUNG (I William sll2000h . 
restricted to one-dimensional problems. It has been extended to 
treat problems in spherical coordinates. Mesh refinement is imple¬ 
mented by bisection of cells, with simple gradient-based refinement 
criteria. It is based on a se cond -order MUSCL scheme using an ex¬ 
act Riemann solver, as in iaile J 199 lh : some improvements have 
been made to the robustness and accuracy of the Riemann solu¬ 
tion in extreme limits. It uses a photon-conservative scheme for the 
int eraction of th e ionizing radiation with the gas, as discussed in 
IWilliamsI 1200 2), but including the obvious modifications for spher¬ 
ical geometry. 

GLIDE is a one-dimensional Lagrangian code, based on the 
same underlying Riemann solver and ionization/radiation solve as 
AQUALUNG (IWilliamsll2002l) . This code uses piecewise-linear re¬ 
construction within mesh cells; however, the interface velocities 
and pressures calculated by the Riemann solver are used to advance 
the zone widths, momenta and energies directly, rather than being 
converted into fluxes. In the calculations presented, the zones have 
been taken to have constant initial spatial width. 

For both of these codes, the hydrodynamical advance was 
made using 7 = 5/3; testing has shown that this makes a mini¬ 
mal difference from 7 close to unity, due the close coupling of the 
radiation step (at both half- and full-step). 


3.2 capreole-C 2 -Ray (G. Mellema, S.J. Arthur) 

The capreole-C 2 -Ray radiation-hydrodynamics code com¬ 
bines the nonrelativistic Roe solver PPM scheme described in 

2 It is not the scope of the present paper to describe each hydrodynamical 
method. We redirect the reader to the cited works for further details. 
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Table 1 . Participating codes in this project. 


Code name 

Code representative(s) 

Kind 

Early Phase 
ID 3D 

Late Phase 

ID 3D 

Mesh motion 

Mesh geometry 1 

AQUILINE 

Williams 

Grid 

X 


X 


Eulerian 

Uniform AMR radial 

capreole-C 2 -Ray 

Arthur 

Grid 


X 


X 

Eulerian 

Uniform cartesian 

FLASH-FERVENT 

Baczynski 

Grid 


X 


X 

Eulerian 

Uniform cartesian 

FLASH-TREERAY 

Wunsch 

Grid 


X 


X 

Eulerian 

Uniform cartesian 

GLIDE 

Williams 

Grid 

X 


X 


Lagrangian 

Uniform radial 

HERACLES 

Tremblin 

Grid 

X 


X 


Eulerian 

Uniform radial 

PION 

Mackey 

Grid 

X 

X 

X 

X 

Eulerian 

Uniform radial (ID) cartesian (3D) 

RAMSES-LAMPRAY 

Haugbqlle & Frostholm 

Grid 


X 



Eulerian 

Uniform cartesian 

RAMSES-RT 

Geen & Rosdahl 

Grid 


X 


X 

Eulerian 

Uniform cartesian 

SEDNA 

Kuiper 

Grid 

X 

X 

X 

X 

Eulerian 

Uniform radial 

SEREN 

Hubber & Bisbas 

SPH 


X 


X 

Lagrangian (SPH) 

Equal mass 

TORUS 

Haworth 

Grid 

X 


X 


Eulerian 

Uniform radial 


1 Refers to the test problem presented in this paper and not to the general capabilities of each code. 


lEulderink & Mellemal (1995) with the radiation transport and pho¬ 
toionization ^ode^ C 2 -Ray (Conservative Causal ray) code de¬ 
veloped by iMellema et al.1 (1200 6b). This code has been used to 
model the expansion of H II regions in turb ulent molecular clouds 
(iMellema et alJl2006al : iMedina et al.ll2014l) and for these applica¬ 
tions the various contributions to the heating and cooling in the 
photoionized an d ne utra l gas are treated by analytical fits, as de¬ 
scribed in lHennev et al.l (2009). For the present benchmark test, a 
simple non-equilibrium prescription for thermal balance has been 
included, which ensures a temperature of 10 4 K in the ionized gas 
and 10 2 K or 10 3 K, as required, in the neutral gas. The calcula¬ 
tions are all performed on a fixed, uniform Cartesian grid in three 
dimensions and limited parallelisation using Open-MP has been 
implemented. 


culated by exploiting the fact that diV is a dimensionless quan¬ 
tity which implicitly takes the dimensionality m of the size of the 
cell Vceii = into account. The ionization rate is then given as 
kion = diVion/ (jih Vceii At). The speed of the I-front expansion 
is captured by enforcing that the change in the abundance of atomic 
hydrogen, \Axn, new — Axu,oid\ = Axh is smaller than 10%. 

Additionally, FERVENT is f u lly c o upled to a r educe d 
chemical network Glover & Clark J 20 I 2 I); I Glover et all d2010h : 
Idover & Mac Low (12007 al ibi): IWalch et al.l d2014l) . It~explicitlv 
tracks the species CO, C + , H + , H and H 2 which allows for a re¬ 
alistic chemical composition and hence a direct way to compare to 
observations. Moreover, it accounts for most heating processes by 
non-ionizing radiation, such as UY photodissociation, vibrational 
pumping and photoelectric heating. 


3.3 FLASH-FERVENT (C. Baczynski) 

FERVENT feaczvnski et al.ll2015l) is a multi-source, Cartesian, ra¬ 
diative transfer code module im plemented in the magnetohydro- 
dynam ical grid code FLASH 4 (IFrvxell et al . 2000 ; Dubev et ahl 
120081) . It is based on an adap tively split ray tracing scheme intro¬ 
duced in I Wise & Abell (2011). Rays are initially cast from a spher¬ 
ically uniform distribution based on the HEALPIX decomposition 
IGorski et al.ll2005l) . They intersect and traverse the domain until 
a splitting criterion is fulfilled, based on the apparent pixel size in 
comparison to the cell face area. In this way a full sampling of the 
mesh domain is guaranteed. Each cell of the mesh is intersected by 
multiple non-aligned rays, which ensures that artifacts from the un¬ 
derlying Cartesian mesh are eliminated. Additionally, the HEALPIX 
sphere is rotated in each timestep to wash out any remaining align¬ 
ment effects. 

It is fully MPI-parallelized by employing asynchronous com¬ 
munication, i.e. as soon as a ray leaves the local domain it is sent 
to its neighbor. This has the advantage that domains without any 
point sources of their own are involved as soon as possible instead 
of having to wait until the neighbor completely sampled its domain. 

Ionising radiation is modeled in a photon conservative, mesh 
resolution independent fashion. Photon conservation is guaran¬ 
teed by an incremental attenuation of the photon flux diV = 
N (1 — e -dr ), where the flux N that enters a neighboring cell 
i + 1 is Ni +1 = Ni — d Ni. N is given as an absolute photon 
number N = PAt, with P in s -1 . An ionization rate is cal- 


3.4 FLASH-TREERAY (R. Wtinsch) 

TREERAY (Wunsch et al., in prep.) is a new radiation transport 
method combining a tree code with reverse ray-tracing. I t has been 
implemented into the hydrodynamic code FLASH (IFrvxell et al.1 
2000 , version 4.2.1 used in this work). The TREERAY method has 
several advantages that make it suitable for being coupled with hy¬ 
drodynamic codes. Particularly, (i) it is relatively fast (costs are 
similar to solving for self-gravity), (ii) the calculation time does not 
depend on number of sources (making it ideal for diffuse radiation 
treatment), (iii) ray-tracing has the highest resolution close to the 
point where the intensity is calculated, and (iv) it is relatively easy 
to parallelize this algorithm on distributed memory architectures. 
The disadvantage is that distant sources and regions of absorption 
are smoothed over a larger volume leading to the artificial numeri¬ 
cal diffusion of radiation. 

The algorithm kernel consists of three steps. In the first one, 
sources of radiation and the absorbing gas are mapped onto the oc¬ 
tal tree and the emission and absorption coefficients are calculated 
for each tree node. In the second step, the tree is traversed for each 
grid cell and tree nodes are mapped onto rays going in all direc¬ 
tions. The ray directions are obtained using the HEA LPIX lib rary fo r 
uniform tesselation of the surface of a sphere IGorski et al.ll2005l) . 
and the mapping coefficients are proportianal to the volume of the 
intersection of each tree node with the cone associated with the ray. 
In this way, no absorbing or emitting gas is omitted. In the third 
step, the one-dimensional radiation transport equation is solved 
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along each ray. This computational kernel is called in each hydro- 
dynamic time-step iteratively, until the radiation field converges, in 
order to account for regions irradiated from several different direc¬ 
tions. 

Various physical processes can be implemented into the code 
by providing prescriptions for the absorption and the emission co¬ 
efficients. In this work, we use a simple on-the-spot approximation 
with only one source emitting a constant number of EUV photons 
per unit time and the number of absorbed photons given by the re¬ 
combination rate to higher-than-ground levels. 

3.5 heracles (P. Tremblin) 

heraclesH dGonzalez et al.ir2007l : I Audit et al.ll201lh is a 3E0hy- 
drodynamical code used to simulate astrophysical fluid flows. It 
uses a finite volume method on fixed grids to solve the equations of 
hydrodynamics, MHD, radiative transfer and gravity. This software 
is developed at the Service d’Astrophysique, CEA/Saclay as part of 
the COAST project and is registered under the CeCILL license. 

HERACLES simulates astrophysical fluid flows using a grid 
based Eulerian finite volume Godunov method. It is capable of sim¬ 
ulating pure hydrodynamical flows, magneto-hydrodynamic flows, 
radiation hydrodynamic flows (using either flux limited diffusion 
or the Ml moment method), self-gravitating flows using a Poisson 
solver or all of the above. HERACLES uses cartesian, spherical and 
cylindrical grids. Current ongoing developments include a multi¬ 
grid method and a multi-group scheme for t he radiative transfer. 

The ionization scheme is described in iTremblin et al.l d201 lh 
and was applied to the study of the formation of pillars and globules 
at the interface of H I I regions and turbulent molecular clouds (see 
ITremblin et al.ll2012al ]lbb. 

3.6 pion (J. Mackey) 

PION is an Eulerian magnetohydrodynami cs code tha t solve s either 
the Euler equations of gas dynamics dMackev & LimlEz OlCT) or the 
ideal magnetohydrodynamics equations dMackev & LiirfeOllh on 
a uniform fixed mesh in Cartesian (ID, 2D, 3D), cylindrical (2D 
axisymmetric), or spherical (ID) coordinates. It uses a finit e vol¬ 
ume, shock-capturing integration scheme dFalle et al.l lll998h with 
geometric_source terms to account for curvilinear coordinates when 
needed JFalld 199 1). Abundances of chemical species and tracers 
are passively advected with the flow. PION is written in C++, is de¬ 
signed to be as modular as possible, is parallelised with MPI by 
dividing the domain into N subdomains (N must be a power of 
2), each controlled by one MPI process, and scales well to at least 
1024 cores for 3D simulations. 

Radiative transfer of ionizing and non-ionizing photons from 
point sources and sources at infinity is followed using a short char¬ 
acteristics raytracer to calculate column densities. This is coupled 
to a microphysics module to solve the rate equations for chemi¬ 
cal species and their associated heating or cooling. Scattered and 
re-emitted photons are treated using the “on-the-spot” approxima¬ 
tion, i.e. they are assumed to be re-absorbed locally. A photon- 
conserving formulation of the ionization rates is used and spectral 
hardeni ng of ionizing radiatio n with optical depth is included (fol¬ 
lowing [Meflemaetal] [20063). The non-equilibrium ionization of 
H is integrated together with the thermal evolution of the gas using 

3 http://irfu.cea.fr/Projets/Site_heracles 

4 Here used as ID only 


CVODE dCohen et al.ll996h . a high-order integrator for coupled dif¬ 
ferential equations, set to use backward differencing with Newton 
iteration. The ionization and heating/cooling source terms are inte¬ 
grated explicitly in the finite volume scheme, using an innovative 
algorithm that preserves the second order time-accuracy of the nu¬ 
merical scheme and dramatic ally reduces the computation required 
for a given error tolerance dMackevll2012h . 

3.7 RAMSES-LAMPRAY (T. Haugb0lle, T. Frostholm) 

The radiative transfer module LAMPRAY (Long characteristics 
AMR Parallel RAY tracing) i s implemented into a derivative of th e 
RAMSES cosmological code dTevssie rll2002l : iFromang et al.ll2006h . 
which has been adapted to the detailed study of star formation. 
The code solves the MHD equations on a fully threaded tree (FTT) 
with support for self-gravity. Important additional physics modules 
compared to the original RAMSES code include accreting sink parti¬ 
cles, coupling to the astrochemistry framework KROME, and many 
smaller changes to improve the stability and quality of the HLLD 
MHD solver in the case of supersonic flows. For the D-type test 
the hydrodynamics is solved using a second order MUSCL scheme 
with an HLLD solver and an isotropic monotonized central slope 
limiter. 

The radiative transfer module employs ray tracing directly on 
the adaptive mesh. These are long characteristics rays covering 
the entire computational domain. The technique is made compu¬ 
tationally feasible by, in addition to the FTT, also using a separate 
ray-oriented domain decomposition, in which the radiative transfer 
problem can be solved efficiently in parallel. A photon-conserving 
second-order accurate TSC interpolation between the two domains 
is used to calculate densities and abundances along the ray, and de¬ 
posit the ionization and h eating rates . The ioniz ation solver is based 
on the c 2 -RAY method iMellema et all 12006b). In a given timestep 
first the ray geometry is established and the MPI and sorting keys 
are set up linking effectively cell centers with points along the ray. 
Then the rates in the timestep are found using an iterative method, 
that repeatedly computes the average ionization rate and the corre¬ 
sponding changes in abundance, until the two converge. The result 
is that in principle arbitrarily large time steps can be taken, and the 
method is only limited by the Courant condition imposed by the 
fluid dynamics. 

The ray scheme allows for a completely arbitrary placement 
of rays, and only domains that actually intersect with rays store 
knowledge about them, making the method scalable to thousands of 
cores. The flexible placement of rays in addition makes it possible 
to use a mixture of rays and solve simultaneously for the diffuse 
and point-source components of the radiation field. 

3.8 RAMSES-RT (S. Geen & J. Rosdahl) 

RAMSES -RT0 dRosdahl et al.ll2013l) is a radiation-hydrod ynamical 
exten sion of the RAMSES cosmological code. RAMSES dTevssied 
l2002h solves the equations of gravitational-hydrodynamics with a 
second-order unsplit Godunov solver on an adaptive mesh, using a 
fully thr eaded tre e struc ture. In this paper, the hydro-solver speci¬ 
fied bv lToro et al. l d 19941) has been adopted. The code is fully (MPI) 
parallel. 

In RAMSES-RT, RT is integrated into the structural framework 

5 The RAMSES-RT implementation is publicly available, as a part of the 
RAMSES code (https://bitbucket.org/rteyssie/ramses) 
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of RAMSES and coupled to the hydrodynamics via interactions with 
hydrogen and helium. The RT is solved on the AMR grid with a 
first-order moment method using the Ml closure for the Eddington 
tensor. This strategy has the advantages that the computational load 
is invariant with the number of radiation sources (moment method) 
and the radiation transport solver is local in space (Ml closure), i.e. 
it requires no information outside the local volume when advecting 
the photon fluid between cells. 

The radiation-thermochemistry of hydrogen and helium is 
solved in RAMSES-RT assuming non-equilibrium, where the 
ionised fractions of those species are tracked explicitly as passive 
scalars that are advected with the gas flow. 

For adaptability of the code, and because the implementa¬ 
tion is designed in large part to simulate galactic and extragalactic 
scales, where radiation sources are highly dynamic, the radiation is 
advected with an explicit solver, which is subject to a Courant con¬ 
dition for the time step. The solution to overcoming the extremely 
small implied time-s tep is to use the so-c alled reduced speed of 
light approximation jGnedin & Abell l200Th : the speed of light is 
simply slowed down by some factor, typically two or three or¬ 
ders of magnitude, such that the hydrodynamical time-step length 
is roughly maintained. In the RAMSES-RT tests described here, we 
reduce the speed of light by a factor 10“ 2 compared to the real 
value. 


3.9 SEDNA (R. Kuiper) 


The ionization module named SEDNA is currently under develop¬ 
ment by Rolf Kuiper. The module is coupled to the static g rid ve r- 
sion of the open source MHD code PLUTO jMignone et alJl2007h . 
Up to now, the module works with static, rectangular, uniform and 
non-uniform grids in 1D-3D Cartesian, cylindrical, and spherical 
coordinates. 

Similar to the h ybrid radiation transpor t solver , introduced in 
iKuiper et"akl d2010bh and iKuiper & Klessenl d2013h . the total ion¬ 
izing radiation field is split into a direct ionization component and 
a diffuse radiation field. The ionization by direct irradiation from 
either a point source at the center of the spherical domain or a 
plane parallel flux in Cartesian geometry is computed by a ray¬ 
tracing step along the first coordinate axis. For the secondary radi¬ 
ation field, the module can either make use of the so-called “on- 
the-spot” approximation, or the ionization of EUV photons cre¬ 
ated by direct recombination into the ground state can be computed 
using the flux-limited diffusion approximation. Ongoing develop¬ 
ment of the module i s bas ed on the numerica l code descrip tions by 
lYorke & Wel3 i ll 9961) and lRichline & Yorkd 1 199112000) as well 


Kudritzki et al.1 ( 19881) . The i 


as the lecture course material in 
ization module SEDNA can be coupled to the radiation transport 
solver mentioned above to allow the determination of gas and dust 
temperature in dusty ionized and non-ionized regions, respectively. 
Future development might include dust scattering and FUY radia¬ 
tion fields, which contribute to the photo-heating of gas as well as 
to Carbon ionization. 

The main purpose of the ionization module SEDNA will be 
its application in modeling the formation and feedback of mas¬ 
sive (proto-)stars. In this sense, the module comprises an exten¬ 
sion to previous n u merical stud i es regarding radiation pressure 
(IKuiper et al.1 l2Q10l 1201 ll I20H IKuiper & Yorkel l2013ah . stellar 
evolution ( Kuiper & YorkeluOl 3bl) . and protostellar outflow feed¬ 
back (IKuiper et al.l2015l) 


3.10 seren (D. Hubber, T.G. Bisbas) 

SEREN (iHubber et alJl201lh is an SPH code designed for star for¬ 
mation, planet formation and star cluster problems. SEREN uses a 
conservative ‘grad-h’ SPH implementation to model the hydrody¬ 
namics, a Barnes-Hut tree to model self-gravity, sink particles to 
model accretion onto protostars and severa l algorith ms to model 


prot ostellar and stellar f eedback (i.e. lStamatellos et al.1 

Bisbas et al . ( 2009 ) used the HEALPIX algorithm (iGorski et al.1 
120051) to tessellate all surrounding directions into discrete vectors of 
equal solid angle that cover the whole surface. In that method, ion¬ 
ising radiation is propagated along each HEALPIX ray, calculating 
the position of the ionisation front and assuming ionisation equilib¬ 
rium and neglecting the diffuse field. The trapezium method is used 
to calculate the density at various ‘evaluation points’ along the ray 
following the radiation. The distance to the next evaluation point is 
given by A r = fih where f± is a dimensionless constant of or¬ 
der unity and h is the smoothing length calculated at the previous 
evaluation point. At the final evaluation point, a bisection method 
is used to accurately determine the location of the final ionisation 
front. The temperature of particles is smoothed around the ionisa¬ 
tion front to prevent the two fluids (hot_and cold) from becoming 
separated and forming a gap (e.g. see |Pricell2008l) . The angular res¬ 
olution of the ionising radiation can be refined at any point by split¬ 
ting a ray into four child rays, which then compute the rest of the 
ionisation integral independently. A ray is split when the width of 
the ray cone exceeds some fraction of the local smoothing length. 
This matches the ray resolution to the gas particle resolution. 

SEREN now contains an updated version of this algorithm with 
two main improvements and optimisations. 


(i) When creating the HEALPIX tessellation, the particles are di¬ 
vided into linked lists along each ray. When the rays are split to 
improve the resolution, the linked lists are also split along each 
ray amongst the four child rays. When the radiation is propagated 
along each ray, the step-size between evaluation points is taken as 
the minimum of the smoothing lengths between the previous and 
next particles. This ensures that the next evaluation point does not 
‘over-shoot’ if there is suddenly a dense region such as a high den¬ 
sity clump or a shock. This means the bisection iteration to find 
the ionisation front is now performed more accurately with fewer 
steps. 

(ii) When walking along the linked lists, if the previous and next 
particles have very similar properties to within some tolerance, then 
the density can be extrapolated from the particles rather than per¬ 
forming another (expensive) tree-walk. This speeds up the calcula¬ 
tion, depending on the chosen tolerance. 


Both of these optimisations lead to faster run times and a more 
accurate determination of the location of the ionisation front. For 
roughly uniform density distributions, both methods give the same 
results. For strong density contrasts, particularly near the ionisation 
front (which can be common later in the simulation), the results 
may diverge due to the different accuracies of finding the ionisa¬ 
tion front position. We perform all ionisation tests using the latter 
implementation only. 


3.11 TORUS (T. Haworth) 

TORUS is principally an Adaptive Mesh Refinement (AMR) Monte 
Carlo radiation transport code capable of continuum, atomic line, 
non-FTE molecular l ine radiation transp ort (e.g. lHarrie ] '2000 : 
iKurosawa et all[2006)1 : iRundle et al.ll201(l) and most recently cou- 
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pled radiation transport and photodissociation region chemistry 
(Bisbas et al. 2015, submitted). It was developed to treat radiation 
hydrodynamics pr oblems using operator spli t photoionization and 
hydrodynamics bv lHaworth & Harries ( 2012 ) as follows. 

TORUS uses a flux conserving, finite difference, total variation 
diminishing hydrodynamics scheme. It uses the Superbee flux lim¬ 
iter (iRoell 19851) and includes a Rhie-Chow interpo lation scheme to 
prevent odd-even decoupling (iRhie & Chowll 19831) . 

Monte Carlo photoionization solves for the time-averaged en¬ 
ergy density dU in each cell by propagating constant-energy e 
packets of photons over the computational grid and counting their 
path lengths l through each cell of volume V 


dU = 


47 tJ u 
c 


dv — 


e 1 
~d\tV 


& 


( 20 ) 


where e is the total source luminosity divided by the number of 
photon packets used. 

The estimated energ y dens ities are used in the photoionization 
equilibrium equation (lOsterbrockll 19891) to solve for the ionization 
balance, which in terms of Monte Carlo estimators is 

n(X i+1 ) = € v la^C) 

n(X l ) A tVa(X l )n e ^ hv 

where n(X' 1 ), a(X l ), a u (X l ) and n e are the number density, re¬ 
combination coefficient and absorption cross section of ion X 1 and 
the electron density respectively. The ionization fraction yields the 
temperature and therefore pressure under the simple thermal pre¬ 
scription used in this paper. To remain consistent with the other 
codes the OTS approximation is employed by terminating the prop¬ 
agation of a photon packet after its first absorption and by using the 
case B recombination coefficient. Using operator splitting, hydro¬ 
dynamics and photoionization calculations are performed sequen¬ 
tially until the simulation end time. The advantage of this approach 
is that many complex features usually only treated by dedicated ra¬ 
diation transport codes c an be included in RHP ap plications (e.g. 
the diffuse radiation field jHaworth & Harriesl2012f) . The disadvan¬ 
tage is that it is computationally expensive, but this is overcome 
owin g to the e fficient parallelization of Monte Carlo radiation trans¬ 
port ( Harriesll2015l) . 


4 THE D-TYPE BENCHMARK TEST 

4.1 Initial Physical conditions 

For the purposes of the test we will use a simplified isothermal 
equation of state for both the ionized and the neutral medium. The 
treatment of the interface where the two media meet is left to each 
numerical method. The test is purely hydrodynamical i.e. it in¬ 
cludes no gravity or magnetic fields. Due to the nature of Eqn. 
we will run two tests: i) to examine the early phase of the D-type 
expansion i.e. where the Spitzer approximation (Eqn. |9j is applica¬ 
ble and ii) to examine the later phase of the D-type expansion, i.e. 
where the Raga expression (Eqn. n is applicable. 

In both early and late phase tests, we consider a uniform re¬ 
gion containing pure hydrogen and an ionizing source emitting 
A/’ LyC = 10 49 photons per second placed at the origin. The po¬ 
sition of the source defines the origin of a Cartesian coordinate sys¬ 
tem. The density of the gas (hydrogen only) is taken to be p 0 = 
5.21 x 10 21 g cm 3 . The temperature and the sound speed of the 
ionized gas are taken to be 7] = 10 4 K and a = 12.85 km s -1 
respectively. 


4.1.1 Early phase 

For the early phase expansion test the temperature of the neutral 
gas is initialized to T 0 = 10 2 K corresponding to a sound speed 
of c 0 = 0.91 km s -1 . For the SPH calculations, we assume the 
dense gas takes the form of a cloud of finite radius, with total mass 
of M c i = 640 M©. The radius of the cloud is therefore R c \ = 
1.257 pc, or R c \ = 4Rst- We evolve the simulation until t en a = 
0.141 Myr. At this time the ionization front has reached the edge 
of the domain (or the cloud in the SPH runs). The upper panel of 
Fig. 0 plots the analytical equations described in Section [2] for the 
early phase expansion. 


4.1.2 Late phase 

Equation [14] gives the stagnation radius as a function of the initial 
Stromgren radius. 

/ . \ 4/3 

-^STAG — ( ~ J Rst 5 (22) 

As we discuss in T4.2.1I the number of SPH particles required for 
the late phase test may lead to prohibitively high computational 
expenses (see Eqn.[27]in particular where there is a dependency on 

(f 1 ) )• We will therefore adopt a ‘neutral’ temperature of T 0 = 

10 3 K while keeping p Q = 1. The corresponding sound speed is 
then c 0 = 2.87kms _1 . As for the early-phase calculations, the 
initial Stromgren radius is Rst = 0.314 pc. The stagnation radius 
K ST ag — 2.31 pc is obtained at t STAG ^3.0 Myr which defines 
the time we terminate the simulation. 

For the late phase SPH calculations, the cloud has mass M c \ — 
8 x 10 3 M© and R c \ = 2.91 pc. The lower panel of Figure Q]plots 
the analytical equations described in Section [2] for the late phase 
expansion. 


4.2 Numerical setup and configuration 

4.2.1 SPH setup 

We choose a slightly bigger cloud radius than the stagnation ra¬ 
dius to avoid the effect of an expanding shell in vacuum when 
t = t STAG . Let g be this additional factor determining the extra 
size of the cloud radius, 


i? c i = gR STAG (23) 

and let f = 9 {%) • 

Suppose that R c \ = fR st , where R c \ is the radius of the 
spherical cloud, R st is the Stromgren radius and / > 1 is a user- 
defined factor. / describes the size of the Stromgren sphere in com¬ 
parison with the size of the spherical cloud. 

The density at t = 0 is homogeneous throughout the cloud. 
Therefore, 


/ 3 4c" 2 

M c \ a B 


(24) 


Let iVgp H be the number of SPH particles consisting the 
Stromgren sphere. All SPH particles have the same mass, ra SPH . 
Then 


3-^St _ ^SPH m SPH 

47Ti?3 t “ 471 -Rl 


( 25 ) 
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Figure 1. Analytical expressions of Eqns. E (Raga-I), E] (Raga-II), 0 
(Spitzer), and 1131 ( Hosokawa-Inutsukal for the ‘early phase’ (upper panel) 
and for the ‘late phase’ tests. The x-axis is the time (Myr) and the y-axis 
the position of fronts (pc). 


and 


po — 


3M c i 


3 n: 


SPH '' "SPH 

4 ' 


(26) 


Combining the above, we obtain = f 3 N gp H . This equa¬ 

tion gives the total number of SPH particles for a specific size and 
resolution of the Stromgren sphere. In all cases we will use a pre¬ 
settled particle distribution i.e. a ‘glass’ structure, which minimizes 
(but does not eliminate) the numerical noise. 

For the test examining the early phase we will use f — 4 and 
iV s S pH = 10 4 . Therefore the total number of SPH particles is taken 
to be iVgp H = 6.4 x 10 5 particles. We also use ra SPH = 10 -3 M©, 
which therefore implies M c \ = 640 M©. 

For the test examining the late phase, the total number of SPH 
particles is: 

<h=S 3 (J)\ S p„- (27) 


Using g = 1.26 (because at this value g 3 ~ 2) and N^ pn — 10 4 
we obtain a total of N^ l pH = 8 x 10 6 SPH particles, each one 
having mass ra SPH = 10 -3 M©. 


Table 2. Number of cells used by the ID participating codes for both early 
and late phase tests (Test 1-ID and Test 2-ID respectively). 


Code name 

No. grid cells 

AQUILINE 

32768 

GLIDE 

378 

HERACLES 

16834 

PION 

131072 

SEDNA 

16834 

TORUS 

1024 


4.2.2 3D grid setup 

The radiation source is placed at the origin, as before. The three- 
dimensional Eulerian calculations are run for a single octant with 
a spatial resolution of 128 3 grid zones. Reflecting boundary condi¬ 
tions are used in the negative directions, and zero gradient bound¬ 
ary conditions in the positive directions. An octant runs from 
{x, y, z} G [0, 3.874 x 10 18 ] cm, corresponding to 4Rst (Rst = 
0.9685 x 10 18 cm = 0.314 pc). 

For the late phase test the simulation was set up in a 
similar way, but with the simulation domain now extending to 
1.26i? S TAG = 2.91 pc to agree with the SPH setup (Sec. 14.2.lT ). 
The ISM density is the same, but the neutral gas temperature is ini¬ 
tialized to 10 3 K instead of 10 2 K, so the pressure ratio between 
ionised and neutral gas is reduced to 20 in the initial conditions. 


4.2.3 ID simulations 

For ID simulations, contributors were encouraged to use any al¬ 
gorithm available for solving spherically symmetric problems. The 
simulation domain was typically set to a somewhat larger value 
than 4Rst so that the shock front remained inside the grid until the 
finishing time of the simulation. For the early phase test, a maxi¬ 
mum radius of 1.5 pc was sufficient, while for the late phase test the 
maximum radius was taken to be 5 pc. Numerical schemes used in¬ 
clude a uniform grid (PlON,SEDNA, HERACLES, TORUS), adaptive 
Eulerian mesh (AQUILINE), and a Lagrangian mesh (GLIDE). The 
different schemes resulted in a significant range of spatial resolu¬ 
tions from code to code. 


5 RESULTS 

5.1 ID runs 

5.1.1 Early phase test (Test 1-ID) 

The six ID hydrodynamical codes use a variety of grids specified 
in TableQ] In these runs, a variety of resolutions have been studied. 
Where practicable, runs were performed at increasing resolution 
until the results converged and the resolution requirements turned 
out to be different for different codes. In this section, we show only 
the highest resolution results, in order to provide a validated refer¬ 
ence solution to which other results may be compared (see Table |3 
for the resolution each ID code used). 

Figure |2] shows the density as a function radius for the six ID 
codes participating in the Test 1-ID. In this figure we plot four dif¬ 
ferent evolutionary times, i.e. t = 0.005, t — 0.02, t — 0.08 and 
t = 0.14 Myr. We also indicate with vertical lines the position of 
the extent of the H II regions as predicted by the various approx¬ 
imations presented in Section |3 All codes agree very well with 
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Figure 2. Density plots of the ID codes for the early phase test. The panels correspond to different evolutionary times. The a>axis in these panels is the 
distance from the ionizing source (pc) and the y-axis is the density (gem -3 ). We plot snapshots at t = 0.005Myr (top left), t = 0.02Myr (top right), 
t = 0.08 Myr (middle left), and t = 0.14 Myr (middle right). We find very good agreement between all codes. We see that at t = 0.005 Myr the codes 
reproduce the position predicted by the Spitzer equation (Eqn.[9j, whereas at the end of the simulation, t = 0.14 Myr, the maximum density is located further 
than the Hosokawa-Inutsuka equation (Ean.fT3V Further details are discussed in £15.1.11 


each other; the peak density and shell width vary between codes 
because of the different spatial resolutions used. At early times the 
position of the thin shell agrees with the Spitzer solution (Eqn. [9} 
but then moves towards the Hosokawa-Inutsuka (Eqn. solution 
even overtaking it at late times. At t = 0.08 and 0.14 Myr it is clear 
that Eqn.[9]no longer correctly describes the evolution. At this stage 
all codes at least partially resolve the shell width. 

The two panels in Fig. [5] show the position of the the ioniza¬ 
tion front, defined as the position at which x\ — 0.5, compared to 
the Spitzer (Eqn. [9} and the Hosokawa-Inutsuka (Eqn. QA} approx¬ 
imations. The lower panel shows the relative error of each code as 
compared to those two analytical approximations. This relative er¬ 
ror is ~ 8% for the Spitzer and < 1% for the Hosokawa-Inutsuka 
approximation at t > 0.08 Myr, and so we conclude that all par¬ 
ticipating codes agree with the Hosokawa-Inutsuka approximation 
and none with the Spitzer formulation at this level of accuracy. 

In Table [3] we show (columns 2 and 4 respectively) the mean 
position of the ionization front (.Rif) at several times and the stan¬ 
dard deviation, a, as obtained by the six different ID codes. We 
find that the position r is remarkably similar between the codes, as 
the relative difference is < 1% in all cases, and that this difference 


decreases with time. In the same table we show (in column 6) the 
mean maximum density reached, (p ma x), as well as the respective 
standard deviation (in column 7). As before, better agreement is 
obtained as t increases. 

At later times, the H II region expansion slows down monoton- 
ically, and so in reality the shock Mach number (hence compression 
factor and maximum density) also decreases monotonically. The 
shell mass and thickness also increase monotonically with time as 
more mass is swept up. This means that the quoted (p ma x) is cer¬ 
tainly an underestimate compared to the analytic solution at early 
times, and gradually reaches the correct value at later times. This 
can be seen in the first panel of Fig.|2] where all of the codes have 
different peak densities because of the differing spatial resolution. 
The largest density achieved (by GLIDE) is significantly larger than 
the mean value quoted in Table [3] 

5.1.2 Late phase test (Test 2-ID) 

Figure|4| shows the density as a function of radius for the six partic¬ 
ipating ID codes in the Test 2-ID for different times, i.e. t = 0.05, 
t = 0.2, t = 0.8 and t — 3.0 Myr. The vertical lines correspond 
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Table 3. This table shows results from the ‘early phase’ test (see I5XTV 
Column 1 shows the time. Columns 2 and 4 show the mean position of the 
ionization front, (-Rif) with the corresponding standard deviation, a, for 
ID codes. Columns 3 and 5 are the respective for 3D codes. Columns 6 
and 7 show the mean maximum density of the shock front (p ma x) and its 
standard deviation, a, for ID codes only. 


t 

(Rif) 

cr/10 

-3 

(p max )/10 19 

cr/ 10- 19 

(Myr) 

(pc) 

(PC) 


(g cm -3 ) 

(gem -3 ) 


ID 

3D 

ID 

3D 

ID 

ID 

0.005 

0.373 

0.378 

2.059 

2.76 

5.7 

3.6 

0.01 

0.430 

0.430 

3.180 

2.96 

5.9 

3.0 

0.02 

0.536 

0.530 

2.488 

3.26 

6.1 

3.3 

0.04 

0.717 

0.707 

3.848 

3.59 

4.2 

1.4 

0.08 

1.009 

0.996 

3.679 

4.35 

2.6 

0.4 

0.14 

1.343 

1.324 

2.848 

4.29 

1.7 

0.2 


Early phase; ID codes (Test 1-ID) 
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Early phase; ID codes (Test 1-ID) 



Figure 3. Position at which the ionization fraction is X[ = 0.5 as calculated 
by the ID participating codes. The upper panel compares the codes with 
the Spitzer (Eqn. [9j solid gray line) and the Hosokawa-Inutsuka (Eqn. m 
dashed gray line) equations. The lower panel shows the realtive error found 
between the numerical solution and either of the two analytical equations 
(solid lines correspond to the errors due to the comparison with the Spitzer 
equation, while dashed lines due to the Hosokawa-Inutsuka equation). 


Late phase; ID codes (Test 2-ID) 
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Figure 5. Comparison between Eqns. [8] (Raga-I), and fill (Raga-II) for the 
ID participating codes. The black dashed line corresponds to the Star- 
Bench Eqn. [28] We additionally plot the Spitzer and the Hosokawa- 
Inutsuka (HI) expansion laws for reference. 


to the positions predicted by the analytical expressions discussed in 
Section [2] Table [4] gives the mean position of the ionization front 
with the standard deviation (columns 2 and 4) as well as the mean 
and standard deviation for the maximum density reached (columns 
6 and 7). The relative difference around (Rif) is < 2% in all cases. 

In this late phase test, the position of the ionization front stag¬ 
nates as predicted by Eqn. [8] (Raga-I) since the H II region is 
in pressure equilibrium with the pressure of the neutral medium. 
Indeed, as shown in Table 0] the mean stagnation radius (mea¬ 
sured at t = 3.0 Myr) as given by all six participating codes is 
{Rif} = 2.359 pc (standard deviation a — 0.032 pc or ~ 1.3% 
from this mean value) which gives a relative error of ~ 2.2% when 
compared to Eqn. [14] 

The shell thickness is a function of resolution, which is why 
the thickness obtained by TORUS is larger than the other results 
(i.e. 1024 cells, compared to a maximum of 40960 used in PION). 
The shell develops rapidly and becomes geometrically thick by 
t — 0.05 Myr. As seen in Fig. [4] the thickness is remarkably sim¬ 
ilar as calculated by all six participating codes. The shell becomes 
thicker over time as observed for t ^ 0.2 Myr. The leading shock 
eventually becomes a detached expanding compression wave, and 
the value of the maximum density p ma x drops down towards ~ p 0 
which is the density of the undisturbed neutral medium (i.e. at the 
beginning of each simulation). 

Figure [5] shows the position of the ionization front of all codes 
compared with Eqns[8] (Raga-I) and \TT\ (Raga-II). For reference, 
we also plot the Spitzer and Hosokawa-Inutsuka equations. As first 
shown by iRaga et all d2012bl) . the codes initially follow the posi¬ 
tion predicted by Eqn. \TT\ however they ‘relax’ towards the stag¬ 
nation radius expected from pressure balance predicted by Eqn. d 
The position of the ionization front temporarily expands beyond the 
stagnati on radius before relaxing inwards. According to [Raga et al.l 
(I2012bh . this ‘overshoot’ is due to the inertia of the H II region. 

Here, we provide a semi-empirical equation (which we call the 
‘StarBench equation’) predicting the position of the ionization 
front for this late phase test for t G [0, 3] Myr. The StarBench 
equation has the form: 


Rsb — R11 + /sb(R-i — R^ii), 


(28) 
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r [pc] r [pc] 


Figure 4. As in Fig. [2] for the late phase test. We plot snapshots at t = 0.05 Myr (top left), t = 0.2 Myr (top right), t = 0.8 Myr (middle left), and 
t = 3.0 Myr (middle right). 


where 1Z\ and IZu are the numerical solutions of Eqns[8]and[TT]re- 
spectively, and /sb is assumed to be a time-dependent dimension¬ 
less factor. For the purposes of the present paper we approximate 
fsB with the following form 


fsB = 1 - A exp 


t 

Myr 


(29) 


where A = 0.733. This equation is plotted in Fig. [5] and [7] and 
shows very good agreement with all ID and 3D simulations. We 
propose that coders who benchmark their algorithm against this late 
phase test should use Eqns[8] \TT\ [28] and [29] with the indicated 
constant for t > 0.05 Myr. 


5.1.3 Summary of ID results 

Overall we find good agreement between all six ID participating 
codes. In order to achieve this we require such a high resolution 
that a 3D calculation would be prohibitively expens ive C$> 10 4 
cells in each dimension for a uniform grid, see also iBisbas et"aD 
120091 . for the relevant discussion in SPH). Although all codes are in 
excellent agreement at sufficiently fine mesh scales, we find limited 
agreement with the existing and widely used analytic forms. We 
attribute this to the approximations made in deriving these analytic 
forms. Many of these approximations are only valid when applied 


Table 4. As in Table [3]for the ‘late phase’ test. 


t 

{Rw) 

cr/10 

-3 

( Pmax )/10 20 

cr/ 10- 21 

(Myr) 

(pc) 

(pc) 

(gem- 3 ) 

(gcm“ 3 ) 


ID 

3D 

ID 

3D 

ID 

ID 

0.05 

0.771 

0.733 

8.685 

2.606 

3.90 

4.0 

0.1 

1.063 

1.025 

4.356 

2.741 

2.43 

1.8 

0.2 

1.460 

1.420 

6.772 

2.645 

1.57 

0.5 

0.4 

1.934 

- 

6.366 

- 

1.05 

0.4 

0.8 

2.375 

2.349 

7.450 

6.739 

0.73 

0.8 

1.6 

2.515 

2.506 

10.07 

5.854 

0.58 

0.6 

3.0 

2.359 

2.424 

31.71 

1.534 

0.53 

0.3 


to the early phase thin shell expansion. We provide a tuned fitting 
formula for which all codes mutually agree to < 0.5% error at all 
times. 
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Early phase; 3D codes (Test 1-3D) 




Figure 6. As in Fig.[3]for the 3D participating codes. 


5.2 3D runs 

5.2 .1 Early phase test (Test 1 -3D) 

In Appendix lAl we present (Figs. lA4l and IA51) snapshots of slices 
through the density distributions at various times during the evo¬ 
lution of the H II region for the eight 3D participating codes. In 
general, the agreement in structure and front position between the 
widely different codes is excellent. 

There are, however, some significant differences. At t — 
0.005 Myr, the peak density of the shell as simulated by the SPH 
code SEREN is significantly lower than the other codes and the shell 
is broader. However, the agreement improves at later times. If the 
number of SPH particles is increased, this leads to a better agree¬ 
ment suggesting that a higher particle count is required to achieve 
comparable resolution to the grid-based codes. The number of SPH 
particles used was tuned to match the overall expansion of the H II 
region (i.e. 10 4 particles for the initial Stromgren sphere) rather 
than to resolve the shell structure. However, the total particle count 
(6.4 x 10 5 ) is comparable to the number of zones in the Eulerian 
calculations (128 3 ). 

As time progresses, evidence of instability appears in some 
of the calculations. In particular, the instabilities are particularly 
prominent in FLASH-TREERAY and CAPREOLE-C 2 -Ray, which 
both use the PPM algorithm to steepen discontinuities. It should 
be mentioned here that in these particular cases the instabilities 



t [Myr] 


Figure 7. As in Fig.[5]for the 3D participating codes. 

develop first in the places where the shock in the neutral gas is 
propagating parallel to the grid and that if a different cut through 
the 3D datacube is taken, the instabilities are not in evidence. This 
strongly suggests that for these t wo codes the instability is of the 
odd-even type, first described by 1 Quirk (11994 1. These stabilities 
can be cured by introducing extra numerical diffusion into the PPM 
scheme when a strong shock is detected parallel to the grid direc¬ 
tion, or using a hybrid Riemann solver, which switches from the 
usual Riemann solver (e.g. the Roe solver used by CAPREOLE-C 2 - 
Ray) to a more diffusive Riemann solver such as an HLLE solver, 
inside shocks. 

Figure [6] shows the comparison of the analytical equations, as 
we discussed above, to the eight participating 3D codes. Here, we 
find that all codes also follow the expansion law as expressed by 
the Hosokawa-Inutsuka Eqn.[l3] with error < 4% at t > 0.08 Myr 
in contrast to > 5% when compared to the Spitzer Eqn.[9] Table [3] 
shows the mean position, (Rif) (column 3), of the ionization front 
and the standard deviation (column 5) for all these codes which we 
find to agree to within 3 — 7% around (Rif). As noted earlier, all 
codes agree with the Spitzer approximation only for t < 0.07 Myr, 
however based on the ID simulations, we argue that this value of 
time will be decreased down to ~ 0.01 Myr by increasing the res¬ 
olution substantially. 


5.2.2 Late phase test (Test 2-3D) 

Figures lA4l and lA5l show cross-section density plots at times t = 
0.05, t = 0.2, t = 0.8, and t = 3.0 Myr. We find that all codes 
are in excellent qualitative agreement until t ~ 0.8 Myr. From 
t « 0.8 Myr, most of the co des s how some level of D-type insta¬ 
bility similar to that shown in I Williams (2002 ). This has the largest 
magnitude in PION, for which the resulting Reynolds stresses ap¬ 
pear to be causing the ionization front to expand to a larger radius 
than found in the other codes. 

In the case of the SPH code SEREN, the shell becomes prone 
to the tensile instability lMonaghanll2000l) and it creates high den¬ 
sity contrast fluctuations in its interior. Eventually at t = 3.0 Myr 
the shell has been completely detached and the hot ionized medium 
is bounded by vacuum. One would expect that this ionized region 
would expand rapidly into the external vacuum as long as the shell 
is detached, however it appears that because there are no SPH parti- 
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cles in this region the SPH summations do not sample it and hence 
the edge of the region functions as a smooth flow boundary condi¬ 
tion. This may be a similar issue to the known co ntact discontinuity 
problem in SPH (lAgertz et al.ll2007l : IPricdl2008h . 

Figure [7] shows a comparison of the position of the ionization 
front as calculated by all participating 3D codes and the analytical 
solutions. We additionaly plot the StarBench Eqn.[28] As we ex¬ 
plained in 115.2.11 none of the codes reproduce any of the analytical 
solutions. Table |4] shows the mean values of the distance of the ion¬ 
ization front and the corresponding standard deviation. All codes 
are in very good mutual agreement within <5%. 


6 SUMMARY & CONCLUSIONS 

We have studied the standard radiation hydrodynamical test case of 
an H II region expanding in an initially uniform density medium. 
Twelve distinct codes participated in this StarBench comparison 
test. We examined two scenarios of the initial ‘early phase’ expan¬ 
sion and the Tate phase’ relaxation to pressure equilibrium with the 
external medium. 

The early phase test shows that the Hosokawa-Inutsuka ap¬ 
proximation (Eqn. d which results directly from the equation of 
motion of the expanding shell and which is in fact a second order 
ODE, agrees with the numerical results. In contrast, the Spitzer ap¬ 
proximation (Eqn. [9j which results from the assumption that the 
thermal pressure of the H II region matches approximately with the 
ram pressure, and which is a first order ODE, underestimates the 
numerical results by a small-but-significant factor (~ 8%), because 
it does not include the effect of the inertia of the material entrained 
in the shell. 

In the late phase test, we tested the codes against Eqns(8] 
(Raga-I) andQj](Raga-II) which include the pressure of the undis¬ 
turbed medium acting on the H II region and are generalizations of 
Eqns[9] (Spitzer) and |T3] (Hosokawa-Inutsuka) respectively. With 
time, this pressure becomes approximately equal to the thermal 
pressure within the H II region, at which point the system reaches 
equilibrium and the expansion of the ionization front is halted. Our 
benchmarking test showed that all participating codes start by fol¬ 
lowing the expansion law as expressed by Eqn.Qj](Raga-II), while 
at later times the H II region stagnates at the position predicted by 
static pressure balance. 

Since the codes do not agree with any of the analytic ex¬ 
pressions, we have developed an analytic parametrization (which 
we call the ‘StarBench equation’) that describes the early and late 
phase expansion of an H II region to within < 2% at all times which 
we recommend be used for future code validations, exercises and 
analytical studies of massive star-forming regions. 

The structure of the expanding H II region is overall in very 
good agreement between all the contributing codes. We have dis¬ 
cussed physics-related issues, such as instabilities, which apply in 
certain case and which may either be physical or numerical. This 
agreement between the codes has improved dramatically as result 
of the comparison test we present in this paper. This comparison 
has increased the confidence of the scientific robustness of the re¬ 
sults of all participating codes. 
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Figure Al. Density plots of the early phase test (Test 1-3D; left) and the 
late phase test (Test 2-3D; right) for all 3D participating codes. The x-axis 
is in pc and the y-axis in p cm -3 . The bars indicate lcr. 


APPENDIX A: FURTHER PLOTS 

Figure lATI shows density plots of the 3D codes for the ‘early phase’ 
test (Test 1-3D; left panel) and the ‘late phase’ test (Test 2-3D; right 
panel). Both panels are similar to the plots of Figs.|2]and|4] The bars 
correspond to lcr standard deviation from all cells/SPH particles. 

Figures |A2]|A5] show cross section (slice) plots of all 3D codes 
for the ‘early’ and ‘late’ phase tests at z = 0 pc. All 3D grid codes 
have used the same spatial resolution (128 cells) and mesh geom¬ 
etry (uniform). For the particular case of the SPH code SEREN, we 
have remapped the SPH partic les on the z — 0 pc slice and have 
assumed a 128 2 resolution (see lPric 320071 for further details). 


This paper has been typeset from a TgX/ LTpX file prepared by the 
author. 
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Early phase; t=0.005 Myr 
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Early phase; t=0.02 Myr 
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Figure A2. Cross section (slice) plots taken at 2: = Opc for the early phase (Test 1-3D) expansion of all participating 3D codes described in 0 x- and y- 
axes are in pc. The logarithmic colour bar shows gas density in g cm -3 . We show snapshots at t = 0.005 Myr and t = 0.02 Myr. 






















































18 


Bisbas T. G. et al. 


Early phase; t=0.08 Myr 
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Early phase; t=0.14 Myr 
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Figure A3. As in Fig. lA2l for t 


0.08 Myr and t = 0.14 Myr. 
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Late phase; t=0.05 Myr 
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Figure A4. Cross section (slice) plots taken at z = 0 pc for the late phase (Test 2-3D) expansion of all participating 3D codes described in 0 x- and y-axes 
are in pc. The logarithmic colour bar shows gas density in g cm -3 . We show snapshots at t — 0.05 Myr and t = 0.2 Myr. 
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Late phase; t=0.8 Myr 
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Late phase; t=3.0 Myr 
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Figure A5. As in Fig. lA4l but for t = 0.8 Myr and t = 3.0 Myr. 




































